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Abstract 

We carry out numerical simulations of vesicle formation based on the density functional 
theory for block copolymer solutions. It is shown by solving the time evolution equations for 
concentrations that a polymer vesicle is spontaneously formed from the homogeneous state. 
The vesicle formation mechanism obtained by our simulation agree with the results of other 
simulations based on the particle models as well as experiments. By changing parameters such 
as the volume fraction of polymers or the Flory-Huggins interaction parameter between the 
hydrophobic subchains and solvents, we can obtain the spherical micelles, cylindrical micelles 
or bilayer structures, too. We also show that the morphological transition dynamics of the 
micellar structures can be reproduced by controlling the Flory-Huggins interaction parameter. 

1 Introduction 

Amphiphilic block copolymers, which consists of hydrophilic subchains and hydrophobic subchains, 
are known to form micellar structures, such as spherical micelles, cylindrical micelles, and vesicles 
[HE]- To clarify how these structures are self-organised is a basic problem of the kinetics of micellar 
systems. 

Several computer simulations have been available to study formation of a vesicle based on 
the particle models (Brownian dynamics (BD) [3], dissipative particle dynamics (DPD) [4], and 
molecular dynamics (MD) [5]). The spontaneous vesicle formation process from the homogeneous 
state observed in these simulations is as follows: The amphiphilic block copolymers aggregate into 
small spherical micelles rapidly from the homogeneous initial state. The spherical micelles grows 
to larger micelles by collision (cylindrical micelles, open disk- like micelles). The large disk-like 
micelles finally close up and form vesicles. The micelle growth process and the closure process 
are slower than the first spherical micelle formation process. Hereafter we call this process the 
"mechanism I" (see also Figure [lja)). Note that the mechanism I is also supported by Monte Carlo 
simulations [6] or experiments of lipid systems [8] . Strictly speaking, the mechanism obtained 
by the experiments for lipid systems may be different from one for polymer systems. However, 
since the DPD simulations for polymer systems and the BD and MD simulations for lipid systems 
give similar results (mechanism I) , we believe that the vesicle formation mechanism is common for 
polymer systems and lipid systems. We also note that the similar mechanism obtained for polymer 
systems by the experiments of morphological transition of cylindrical micelles to vesicles [9]. 

The field theoretical approach, which has been developed to study mesoscale structures of block 
copolymers [TDJ [TTJ Q~2] EH [HI E3 US] , is considered to be useful for vesicle formation, but most of 
the works have been limited to simulations in thermal equilibrium [T71 [18] . Since these simulations 
ignore realistic kinetics (for example, local mass conservation is not satisfied), the vesicle formation 
process observed in these simulations are different from the mechanism I: The first process is similar 
to the mechanism I. Small spherical micelles are formed rapidly. The spherical micelles then grow 



1 



up to large spherical micelles by the evaporation-condensation like process. The large spherical 
micelles are energetically unfavourable, and thus the large micelles take the solvents into them to 
lower the energy. We call this process the "mechanism II" (see also Figure Qlb)). 

The dynamical simulations for diblock copolymer solutions were carried by Sevink and Zvclin- 
dovsky [T5], using the dynamic self consistent field (SCF) theory. However, in their simulations both 
subchains are hydrophobic and the resulting structures are so-called onion structures [20, 21 , 22lll9j. 
The onion structures obtained by simulations are essentially microphase separation structure in 
the block copolymer rich droplets in the solvents, and qualitatively different from the multilayer 
vesicle structures which is often called as "onion structures" in surfactant solutions. The vesicle 
structures contain solvents inside them and therefore different from these polymer onion structures, 
and one should distinguish the onion formation process from the vesicle formation process. We 
also note that the simulations for onion structures were carried out for weak segregation region and 
this does not agree with previous equilibrium simulations for vesicles, since vesicles are observed 
in rather strong segregation region. (In this work, we use the word "strong segregation region" as 
the region where the minimum and maximum values of the density field is approximately and 
1. One may call such region as the intermediate segregation region.) It should be noted that the 
formation processes of the onion structures by simulations [2Tj fl9] are mainly the lamellar ordering 
process in droplet like regions. Therefore the resulting morphologies (onions) are similar to the 
phase separation patterns in droplets |23j than vesicles. Thus we consider that the onion formation 
process is not the mechanism I. 

Recently, dynamical simulations of vesicle formation have been done by He and Schmid [24], 
using the external potential dynamics (EPD) [25]. They obtained vesicles, but their formation 
process is similar to the mechanism II. Unfortunately the mechanism II is qualitatively different 
from the mechanism I, and this means that their result does not agree with particle simulations 
or experiments (we will discuss the reason for this difference in Section |4~3|) . We expect that the 
mechanisms observed by different simulation methods should be the same. 

In the present work, we apply our previous model, the density functional theory for block 
copolymers [22J[TS], to the dynamics. We carry out numerical simulations for amphiphilic diblock 
copolymer solutions in three dimensions. By using the continuous field model, we show, for the 
first time, that a vesicle is spontaneously formed from a disordered uniform phase. The simulation 
result is consistent with the mechanism I. We also show that we can simulate the spontaneous 
formation process of various micellar structures (such as spherical micelles or cylindrical micelles) 
and the morphological transition dynamics. 

2 Theory 

The dynamics of block copolymer systems are well studied by using the dynamic SCF theory 
[T5] as well as the time-dependent Ginzburg-Landau (TDGL) equation [2B]. The dynamic SCF 
simulations are known to be accurate for from weak segregation region to strong segregation region, 
but they consume memory and need large CPU power. In contrast, the TDGL simulations need 
less CPU power and enables large scale simulations. However, the free energy functionals [10l [11] 
used in the TGDL approach is generally not appropriate for the strong segregation region, since 
the validity of the Ginzburg-Landau (GL) expansion is guaranteed only for the weak segregation 
region where the density fluctuation is sufficiently small [571 El H2] • Actually Maurits and Fraaije 
[2"8] showed that the widely used fourth-order GL expansion model is not sufficient for dynamical 
simulations. To overcome this limitation we need to use free energy functional model such as 
the Flory-Huggins-de Gennes-Lifshitz free energy [29l [30l [31] , which is not expressed in the GL 
expansion form. 

In previous simulations, vesicles are observed for rather strong segregation region in diblock 
copolymer solutions. Thus we can conclude that the use of inappropriate free energy functional 
models qualitatively affects micellar structures in diblock copolymer solutions. We need to use a 
free energy model and a dynamic equation which is valid for the strong segregation region and 
for the macrophase separation. In our previous work [22) we proposed the free energy functional 
which is valid for strong segregation region, that is, valid for vesicles [15] . 
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2.1 Free Energy Functional 

The free energy functional for the system can be expressed as follows by using the soft-colloid 
picture [Ml [S3 131] • 

F = -TS + U = -TS ttaBB -TS CODt + U (1) 

where S and U are the entropy functional and the interaction energy functional, and T is the 
temperature. The entropy functional can be decomposed into two parts; the translational entropy 
functional 5trans and the conformational entropy functional S con f. In the density functional theory, 
U, S'trans and 5 con f are expressed as the functional of the local volume fraction fields. (The local 
volume fraction field is equivalent to the local density field, if the segment volume is set to unity. 
We assume that the segment volume is unity in this work.) 

Each contributions to the free energy functional can be modelled for AB type diblock copolymer 
and solvent mixtures by [22j [T8] 




where ks is the Boltzmann constant, <fii(r) is the local volume fraction of the i segment at position 
r, ipi(r) is the order parameter field (?/>-field) defined as ipi(r) = yj 4>i(r) 31j. The coefficients 
A-ij cLIld Cij are constants determined from the block copolymer architecture (we don't show the 
explicit forms of A+j and Cij here; they can be found in Refs [T8], [22]), /< is the block ratio of the 
block copolymer, b is the Kuhn length, Xij is the Flory-Huggins x parameter, and Q(r — r') is the 
Green function which satisfies [—V 2 + \~ 2 }G(r — r') = S(r — r') [231 [3S], where A is a cutoff length 
and is about the size of microphase separation structures. P(r) is the Lagrange multiplier for the 
incompressible condition (0a(*°) + 0s( r ) + 4>s{ r ) — 1) [36] . 

The conformational entropy ([3]) is expressed in the bilinear form of ■0-field [37l [30l [22] . One 
significant property of the conformational entropy is that it satisfies the following relation. 

Sconf [{<*&(!■)}] = aSconf [{^<(r)}] (5) 

where a is arbitrary positive constant. (It is clear that the conformational entropy in the bilinear 
form of ip satisfies eq ([5]).) We can interpret eq ([5]) as the extensivity of the conformational entropy; 
under the current approximation (using mean field, and neglecting many body correlations), the 
total conformational entropy is sum of the conformational entropy of each polymer chains. Since 
the density profiles of miecellar structures are determined to achieve the low free energy, the 
conformational entropy which does not satisfy eq ([5]) may lead unphysical density profile. We note 
that the Lifshitz entropy [371 130] and the conformational entropy calculated by the SCF theory 
satisfy eq ([5]), while most of previous phenomenological free energy models do not satisfy eq ([5]). 



3 



Substituting eqs © - dl]) into eq (P), we get 



k B T 



" B) J 

f drdr^^WoUMr-r')^)^^') 

i,j (=A,B) 

+ / dr WfAfBC AB ipA (r)i/> B (r) + J2 / dr T l v ^(r)| 2 

i{=A,B,S) 

E / dr^?(r)$(r) + J dr^l [&(r) + tfe(r) + ^(r) - l] 



(0) 



The difference between the free energy functional ([6]) and our previous model [22] is the form of the 
Green function (in the previous theory, the Green function Q{r~r') has no cutoff, — V 2 Q(r — r') = 
S(r — r')). A cutoff for the Green function was first introduced by Tang and Freed |34| . without 
derivation, for the Ohta-Kawasaki type standard GL free energy [TT]. The Ohta-Kawasaki type 
Coulomb type long range interaction term may cause unphysical interaction or correlation in non- 
periodic systems such as block copolymer solutions. The true form of long range interaction can 
be, in principle, obtained by calculating higher order terms in free energy functional (the effect of 
higher order vertex functions) exactly. While we can calculate the higher order terms numerically 
by using the SCF, it is practically impossible to get analytical form of them and take a summation. 

To avoid the unphysical correlation without numerically demanding calculations, in the present 
study, We employ a theory which describes more precisely the properties of the length scale of the 
order of one polymer chain [35]; We modify the Coulomb type interaction by introducing Tang- 
Freed type cutoff [33] , instead of calculating higher order terms exactly. It is reasonable to consider 
that the interaction range (or the cutoff length A) of the long range interaction cannot exceed the 
characteristic size of the polymer chain. Generally it is difficult to determine the cutoff length (one 
way to determine it is to use the SCF calculation or other microscopic calculations such as the 
MD). 

Here we estimate the cutoff length for the simplest case, diblock copolymer melts. For homo- 
geneous ideal state, the cutoff length is considered to be the mean square end to end distance of 
a polymer chain N x ^ 2 b where TV is the polymerization index. On the other hand, for microphase 
separation structures at the strong segregation region, the polymer chains are strongly stretched. 
The periods of the structures are known to be proportional to N 2 / 3 b [TTJ [38] . Thus we have the 
following polymerization index dependence for the cutoff length of diblock copolymer melts. 

\oc{ Nl/2b ( for ideal state ) ( 7 ) 
[ N 2 / 3 b (for strong segregation region) 

We assume that the micellar structures have the similar cutoff length as the melt case. In 
this work we use A ~ N 2 / 3 b for micellar systems. Note that this is a rough estimation and the 
validity is not guaranteed. For example, the proportional coefficient in cq ([7]) is ignored here, or the 
effect of swelling of the hydrophilic subchain in the solution is not taken into account. Therefore 
the validity of this value of A should be tested by comparing with the characteristic scale of the 
resulting phase separation structure. We also note that the dependence of the morphologies to the 
value of A is not so sensitive (see also Appendix [C]) . What important here is that the interaction 
range is not infinite (as the original Ohta-Kawasaki green function) but finite (as the modified 
green function by Tang and Freed). 

2.2 Dynamic Equation 

We employ the stochastic dynamic density functional model [39[ [40], [41] for the time evolution 
equation. The diblock copolymer solutions are expressed as three component systems (hydrophilic 
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subchain A, hydrophobic subchain B, and solvent S). The dynamic equation for multicomponent 
systems can be expressed as follows. 



dt 



Q 0(pi(r) 



(8) 



where 4>i(r) is the concentration of i-th component (i = A, B, S) and £i is a friction coefficient. F 
is the free energy. £,i(r,t) is the thermal noise which satisfies the fluctuation-dissipation relation. 



Ur,t))=0 

^i(r,t)^(r',t')) = —^kBTSijV [&(r)V*(r -r')]5(t-t') 

Si 



(9) 
(10) 



where (...) represents the statistical average and < < 1 is the constant determined from 
the degree of coarse graining (see Ref J41] for detail). One can interpret that the temperature of 
the noise is the effective temperature f3~ 1 T, instead of the real temperature T. We also note that 
£i(r) can be generated easily by the algorithm proposed by van Vlimmeren and Fraaije [42] . It 
should be noted here that the free energy functional in eq ([8]) is generally some kind of effective 
free energy functional, and is not equal to the free energy functional for the equilibrium state [41]. 
In this work we approximate the effective free energy functional as the equilibrium free energy 
functional ([6]). This approximation corresponds to assume that all the polymer chains are fully 
relaxed, and thus this approximation neglects the viscoelasticity associated with the relaxation of 
polymer chains. 

We can interpret the dynamic equation ([8]) as the TDGL equation. In this case, the Onsager 
coefficients (or the mobility) <j)i(r)/Q is proportional to <j>i(r). This is most important for the 
strong segregation region or in the situation that the solute concentration is sufficiently small 
[121 0H [35] . Here we note that the TDGL equation with density dependent mobility and the free 
energy of ideal gases (the translational entropy of ideal gases) reduces to the diffusion equation (or 
the Smoluchowski equation) [45]. This suggests that we should use the Flory-Huggins (or Bragg- 
Williams) type free energy model once we employed the density dependent mobility. Similarly we 
should use the density dependent mobility if we employ the Flory-Huggins type free energy model. 
It is also noted that the TDGL equation with the constant mobility and the Flory-Huggins type 
free energy leads unphysical singular behaviour near <pi{r) = 0. 

For simplicity, we assume that all the segments have the same friction coefficient (Q — C). We 
can rewrite eq ([8]) by using 4>i(r) = ipf{r) as follows. 



dt 



C 



<(r)V 



S(F/k B T) OMr) 



S^(r) dipf{r) 



+ £i(r,t) 



k B T 

2C 
k B T 

2C 



[&(r)V 2 w(r) - /i 2 (r)V 2 ^(r)] + &(r, t) 

^(r)V 2 Mi (r-) - Mr)V 2 Ur) + ^&(r,t) 

k B T 



(11) 



where Hi(r) = 8(F /k B T) / 8tpi{r) is a kind of chemical potential field. It is noted that /J,i{r) 
is not singular at ipi( r ) = and thus we can perform stable simulations (on the other hand, 
S(F I k B T) I S4>i(r) has a singularity and is numerically unstable). 
By introducing rescaled variables defined as 



- k B T 



ti(r,t) 

and substituting eqs (fH?)) and (fHij) into eqs (fTTj) , 
rescaled variables 



2C 
2C 
k B T 



i(r) 



at 



(12) 

6(r,f) (13) 
and (|10[) . we obtain the dynamic equation for 



(14) 



il>i (r ) V 2 ^ (r ) - m (r ) V 2 ^ (r ) + & (r , t) 
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and the fluctuation dissipation relation for the rescaled noise field 

(f<(r,t))=0 (15) 

(&(r, *)&•(/,*')) = -4r x %V • [&(r)W(r - r')] ?) (16) 

It should be noted here that the magnitude of the noise £i(r,i) depends only on The temper- 
ature change corresponds to the change of the time scale (since t cx t/T) and the change of the x 
parameter. However, notice that to return from the to rescaled time t (which is used in the actual 
simulations) to the real time scale t, we have to multiply the factor 2^/ksT. 

We note that we can set ksT/2C, = 1 instead of introducing rescaled variables (t and This 
corresponds to setting the effective diffusion coefficient for the monomer to unity (strictly speaking, 
the half of the effective diffusion coefficient is set to unity) . This change can be done without losing 
generality, as shown above. To return to the real time scale, we have to multiply the factor 2^/ksT 
to the rescaled time. This factor can be estimated from the experimental diffusion data. (We will 
estimate the real time scale in Section l4~2l ) 

It should be noted here that eqs (fl~4|) and (fT6|) implies that the annealing process (the temper- 
ature change process) corresponds to the change of the x parameter, and the magnitude of the 
noise is not changed because there are no other parameters related to T. The magnitude of the 
noise is characterized only by in the rescaled units. 



3 Simulation 

We solve eq (fl~4")l numerically in three dimensions. The chemical potential jUj(r) can be calculated 
from the free energy functional ([6]) whereas thermal noise t) is calculated by the van Vlimmeren 
and Fraaije's algorithm [42 . 

3.1 Numerical Scheme 

Simulations are started from the homogeneous state (<pi(r) — <f>i, where <pi is the spatial average 
of 4>i(rj). Each step of the time evolution in the simulation is as follows; 

1. Calculate the t/'- field ipi(r) and the chemical potential field ^i(r) from the density field <f>i(r). 

2. Generate the thermal noise £i(r,t). 

3. Calculate the Lagrange multiplier P(r). 

4. Evolve the density field <pi(r) by time step At, using ipi(r), //j(r) and £i(r,t). 

5. Return to the stepQ] 

Here we describe the numerical scheme for the simulations in detail. As mentioned above, 
the dynamic density functional equation ([5]) is reduced to the diffusion equation for the Flory- 
Huggins type translational entropy functional. In our free energy (0, we have the Flory-Huggins 
translational entropy term. Thus the dynamic equation (|14p contains the standard diffusion term 
(Laplacian term). 



dt 



2C i V 2 4 (r) + --- (17) 



where we set d = fiCu (i = A,B) and C s = 1. The numerical stability of the diffusion type 
equation can be improved by the implicit schemes. In the previous work [18j we employed the 
alternating direction implicit (ADI) scheme [46j for the Laplacian term. In this work we employ 
the ADI scheme to improve stability, as the numerical scheme for the equilibrium simulations. 

We split the each time evolution steps (time step At) into three substeps (time step At/3). The 
update scheme for the density field from <p\ n \r) = <pi(r,t + nAt/3) to (f>l n+1 ^ (r) is as follows. (For 
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simplicity, here we use the continuum expression for the position r and the differential operators 
d 2 /dx 2 ,d 2 /dy 2 ,d 2 jdz 1 and V 2 . In real simulations we use the standard lattice for the position 
and the center difference operators for the differential operators [46] ) 



2At - d 2 
l -— Ll dx 2 ~ 



1 _ ^Ci — 
3 1 dy 2 



2At~ d 2 

1 ~ Ll ~d^ 



# ) (r)=^° ) (r) + f 



# ) (r)=#(r) + f 



«tfV)=*f°(r) + f 



-2C i ^^(r) + ^(r)V^,;:"' ( rl 



72„(0)/ 



-2ci4 1) W + V'| 1) (r)V 2 / xf ) (r) 



9y 2 



(18) 



(19) 



(20) 



In each steps, the chemical potential fj,\ n \r) and ^ n \r) is regenerated. The numerical Scheme for 
the calculation of the chemical potential (j!- n \r) is found in Ref [IS] • Eqs (fT8|) - (|20|) can be solved 
easily by using the numerical scheme for the cyclic tridiagonal matrix |46| . 

We need to calculate ip^ (r) to calculate the chemical potential /4 (r). Because of the thermal 
noise and numerical error, the condition < (p- n (r) < 1 is not always satisfied. Thus we put 
tpl (r) as 



(0( n) (r)<O) 

r)-'(r) = {l (4 l \r)>l) 

\J 4>[ n \r) (otherwise) 



(21) 



This avoids numerical difficulty associated with negative value of (f>\ n \r). 

The noise field ^i n \r,t) is calculated at each ADI time steps, by using the scheme described 
in Appendix [Al (notice that the size of the time step here is At/ '3 instead of At). 

The Lagrange multiplier P(r) is also updated at each ADI time steps by the following scheme. 
Because of the large thermal noise, the incompressible condition is not always satisfied. Thus we 
use roughly approximated and relatively simple scheme to calculate P(v). If we ignore the terms 
in the chemical potential except for the Lagrange multiplier term, the ADI update scheme can be 
approximately written as 



(22) 



From eq (|22|) we have 



£^ l) (r)VP(r) 



1 



£<^ +1) (r)-£^ n V) 



V 2 P(r) - V 



VP(r) 



(23) 



At/3 

Using the constraint J^i </ > i™ +1 ' ) ( r ) — 1 an d assuming |1 — (r)| -C 1, finally we have 



V 2 P(r) 



1 

JaJ/S) 



(24) 
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Eq ([2"4"| and the long range interaction terms (which contains the Green function Q(r — r')) in 
the chemical potential (J-i(r) can be calculated by several numerical methods for partial differential 
equations. Here we calculate them by the direct method, using FFTW3 [47] . 

3.2 Results 

Parameters for the diblock copolymer are as follows: polymerization index N — 10, block ratio 
Ja = 1/3, fs = 2/3, the spatial average of the volume fraction </> p = 0.2. The cutoff length for 
the long range interaction should be set appropriately. If there are no solvent diblock copolymers 
form strongly segregated microphase separation structure (currently we are interested in the strong 
segregation region). Thus from eq ((7]) we have A oc N 2 / 3 b ~ 4.64. Here we set A = 5. The validity 
of this value is argued later. (See also Appendix ICl) 

Parameters for the solvent are as follows: the spatial average of the volume fraction 4> s — 
1 — 4> p — 0.8. Flory-Huggins \ parameters are xab — 2.5, xas — — 0-5, xbs = 5 (the A monomer 
is hydrophilic and the B monomer is hydrophobic) and the Kuhn length is set to unity (6 = 1). All 
the simulations are carried out in three dimensions. The size of the simulation box is 246 x 246 x 246 
and the number of lattice points is 48 x 48 x 48. We apply the periodic boundary condition. The 
time step is set to At = 0.0025 and the magnitude of noise is = 0.0390625. The simulation is 
carried out up to 2500000 time steps (from t = to t = 6250) and it requires about 20 days on a 
2.8GHz Xeon work station. 

The snapshots of dynamics simulations are shown in Figure [2l The observed vesicle formation 
process is different from the mechanism II, which is observed in the previous thermal equilibrium 
simulations [TTl [18] and EPD simulations [24]. First, spherical micelles are formed (Figure Ufa)). 
Then the micelles aggregate and become large (Figure [5Jb)), and a disk- like micelle (Figure (He)) 
appears. The disk-like micelle is closed (Figure [5Jd),(e)) to form a vesicle (Figure [^f)). This 
vesicle formation process agrees with the mechanism I, which is observed in the particle simulations 

To confirm that our vesicle formation is independent of the random seed, next we perform 
simulations with different random seeds. We peform four simulations here and vesicles are obtained 
in two simulations. Figures and 2] show the snapshots of vesicle formation processes obtained for 
different random seeds (all other parameters are the same as the previous simulation) . The vesicle 
formation processes in Figures [3] and [4] again agree with the mechanism I. While the vesicles are 
not observed for all the simulations, but we consider that the vesicle formation process obtained 
by the previous simulation is not artifact. 

The characteristic size of the phase separation structure should be compared with the cutoff 
length A. The characteristic size of the micellar structure can be found in the two dimensional cross 
section data or the density correlation function. Figure [5](a) and (b) show the two dimensional 
cross section of the final structure (t = 6250, Figure [2](f)). From Figure [5ja) and (b) we can 
observe that the cutoff length A = 5 is comparable to the characteristic size of the micelle bilayer 
structure. Figure [5] (c) shows the radial averaged correlation function of the density field in the 
Fourier space Sab(<}), defined as [15] 



E 



^A(q')/fA-0 B (q')/fB} 2 



SAB{q) = = (25) 

q-A/2<\q'\<q+A/2 

where <f>i{q) is the Fourier transform of the density field and A is the width of the shell in the Fourier 
space (here we set A = 0.5 x 27r). The correlation functions have peaks at g/27r ~ 0.15. Thus the 
characteristic wave length for the corona (the A subchain rich region) and the core (the B subchain 
rich region) is ~ 6.7. While the resolution (the shell width A) is not fine, the characteristic size of 
the structure is comparable with the value of the cutoff length, A = 5. Therefore we consider that 
the value of A used in this simulation is appropriate and valid. (See also Appendix [C| . 

Figure [5] shows various micellar structures obtained from simulations for different volume frac- 
tions (j>p = 0.1, 0.15, 0.25, and 0.3. Simulations are performed up to 1250000 time steps (from t = 
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to t — 6250). All other parameters are the same as those in Figure 03 If the volume fraction of 
the polymer is small (Figure a), (b)) only small micelles (spherical micelles in^a) and disk-like 
micelles in[6jb)) are formed. On the other hand, if the volume fraction of the polymer is large 
(Figure [BJc),(d)) bilayer structures are formed. Vesicles in Figure [5] are observed only for the 
intermediate volume fraction. 

In experiments, the mixed solvent (mixture of organic solvents and water) is often used to 
control the morphology of the micellar structures of block copolymers [55], [SUl EI] - If we assume that 
the mixed solvent is the mixture of common solvent (xas and Xbs is sufficiently small or negative) 
and water, the volume fraction of water mainly affect the interaction between the hydrophobic 
subchain and the solvent. We change xbs to mimic these experiments. 

Figure [7] shows the results of simulations for various xbs {xbs = 0.2, 0.25, 0.3 and 0.35). We can 
observe the spherical micelles are formed for low Xbs and cylindrical micelles are formed for high 
Xbs- This agree with the experimental results qualitatively. Shen and Eisenberg [151 [SO] reported 
that by controlling solvent condition, one can control the morphology of the block copolymer 
micelles. The control of the morphology is especially important for application use such as the 
drug delivery system. Here we mimic the solvent condition change by changing the x parameter 
between the hydrophobic segment and the solvent, xbs- We change the x parameter at t = 6250 
(Figure E^f)), from xbs = 5 to xbs — 2.5 or 3, and perform the simulations up to t = 9375. 
The snapshots at t = 9375 are shown in Figure [8] The vesicle is changed into spherical micelles 
(Figure IH a)) or cylindrical micelles (Figure [UJb)). The resulting morphologies are consistent with 
the morphologies obtained by the simulations started from the homogeneous state (Figure [7](b) 
and Figure COJc)). 

4 Discussion 

From the dynamics simulations based on the density functional theory for amphiphilic diblock 
copolymers, we could observe the time-evolution of spontaneous formation of vesicles. We also 
observed the micellar structures depending on the volume fraction of polymers. The observed 
vesicle formation process is just the same as the results of the previous particle model simulations 
(mechanism I). At the initial stage of time evolution, rapid formation of small spherical micelles is 
observed. After the small micelles are formed, they aggregate each other by collision and become 
larger micelles. Finally the large disk-like micelles are closed and become vesicles. The late stage 
process, collision process and the close-up process, is slower than the initial stage, as the mechanism 
I. 

There are several important differences between our model and the standard TDGL model for 
the weak segregation limit. Here we mention some essential properties of our model. First, our 
model can be applied for the strong segregation region, at least qualitatively. This is especially 
important, since the vesicles are observed in the strong segregation region. 

Second, the interaction between the hydrophilic subchain and solvent should be handled care- 
fully. If the hydrophilic interaction {xas) is too small or both subchains are hydrophobic, the 
system undergoes macrophase separation and separates into block copolymer rich region and sol- 
vent rich region. In such a situation, vesicles are not formed and onion structures [20] [211 E21 [19] 
are formed instead. 

Third, we applied rather large thermal noise to the system. Recently Zhang and Wang JS5] have 
shown that the glass transition temperature and the spinodal line (stability limit of the disordered 
phase) for the microphase separation in diblock copolymer melts are quite close. This implies 
that microphase separated structures in the block copolymer systems are intrinsically glassy at 
the strong segregation region. Without sufficiently large and realistic thermal noise, the system is 
completely trapped at intermediate metastable structures (in most cases, spherical micelles) and 
thus thermodynamically stable structures (vesicles) will never be achieved. As mentioned above, 
the magnitude of the noise is determined by the characteristic time scale and the factor is 
independent of the temperature. Thus the noise is important for mesoscopic systems even if the 
temperature of the system is not high. 
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4.1 Thermal Noise 



Here we discuss about the thermal noise in detail. Intuitively, the large thermal noise is needed to 
overcome the free energy barrier. In the kinetic pathway of usual macrophase separation processes 
(for example, phase separations in homopolymer blends), there are not so large free energy barriers. 
Therefore, in most cases, the simulations for such systems work well with small thermal noise or 
without thermal noise. On the other hand, in the kinetic pathway of vesicle / micelle formation 
processes, there are various large free energy barriers (for example, in the collision and coalescence 
process of spherical micelles). However, it is not clear and established how to determine the 
magnitude of the thermal noise for such systems. 

As mentioned, the thermal noise in eq (jpj) is expressed by using the effective temperature 
/3 _1 T, instead of the real temperature T and /3 _1 is determined from the degree of coarse graining. 
Intuitively this means that the coarse graining procedure decreases the degree of the freedom of 
the field. The similar approach can be found in the van Vlimmeren and Fraaije's theory [42] . They 
employed the noise scaling parameter to express the effective degree of freedom in the simulation 
cell. These two approaches are similar and in fact, we can write the relation between them. 

n = (26) 

where h a is the lattice vector. From eq (|26p we can calculate the noise scaling parameter for our 
simulation. Using the parameters used in Section I3~2l we have ft = (0.5 3 )/0. 0390625 = 3.2. It 
should be noticed that this value of the noise scaling parameter is much smaller than the values 
used in most of previous works (f2 = 100) by Fraaije et al [53l [54j [55] . However, since the thermal 
noise play an important role in the vesicle formation process, we cannot underestimate /3 . One 
may consider that the mean field approximation and the free energy functional ([6]) does not hold 
for low noise scaling parameter systems. It may be true, but as shown by Dean 39], the dynamic 
density functional equation holds even if there are not large number of particles. Thus here we 
believe that the mean field approximation still holds for our systems. 

To show the necessity of the large thermal noise in the simulations, we perform several simula- 
tions with various values of /3 _1 . We expect that the simulations for small boxes are sufficient to 
see the effect of the magnitude of the thermal noise, since the effect is large even for early stage of 
the dynamics. Thus we perform simulations for small systems (system size 126 x 126 x 126, lattice 
points 24 x 24 x 24) with different values of ff- 1 = 0.00390625,0.0390625,0.390625. All the 
other parameters are set to the previous vesicle formation simulations in Section |3~21 (Figure [2]). 
Figure [9] shows the snapshots of the simulations. For small thermal noise case (Figure Ufa), (b)), 
we can observe many micelles with sharp interfaces. For large thermal noise case used in Section 
13.21 (Figure E^c),(d)) we can observe the collision-coalescence type coarsening process of micelles, 
driven by the thermal noise. For larger thermal noise case (Figure G2e),(f)) no structures are ob- 
served. (We consider this is because the thermal noise is too large.) Thus we can conclude that the 
large (but not too large) thermal noise is required. This result is consistent with the proposition 
of van Vlimmeren et al [56] . They carried out two dimensional simulations for dense amphiphilic 
triblock copolymer solutions, with the full (large) noise and the reduced (small) noise and proposed 
that the real systems should have some intermediate noise (smaller than the full noise, but larger 
than the reduced noise). 

It should be noted that simulations with small thermal noise is very difficult to perform for our 
systems, because they are in the strong segregation region. The simulation with the small noise 
shows strong lattice anisotropy. This means that there are very sharp interfaces which often leads 
numerical instability, and thus we should use finer lattices for such systems. On the other hand, 
we do not observe such strong lattice anisotropy in the simulations with large noise. We believe 
that the large thermal noise stirs the density fields and stabilizes the simulation. 

At the end of this section, we note that it is difficult to determine the parameter /3 _1 theoreti- 
cally. The estimation of the magnitude of the thermal noise is the open problem. 
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4.2 Comparison with Experiments 

Next we discuss about the relation between the simulation results and the experiments. We used 
many parameters for the simulations and they should be related to the real experimental parame- 
ters. Unfortunately, because the accuracy of our theory is not high and many approximations are 
involved, it is difficult to calculate the real parameters from our simulation parameters. 

The most important parameters are the Flory-Huggins x parameters. The x parameters are 
experimentally measured or calculated from the solubility parameters. The x parameter between 
the hydrophilic subchain and the hydrophobic subchain is in most cases sufficiently large to cause 
the microphase separation. For example, Bhargava et al [51] calculated the x parameter between 
polystyrene (PS), which can be used as the hydrophobic monomer [2], and water as XPS.water = 6.27. 
Dormidontova [57] measured the x parameter between poly(ethyrene oxide) (PEO) and water and 
reported that the x parameter can be expressed as XPEO.watcr = — 0.0615 + 70/T {T is the absolute 
temperature). Lam and Goldbeck-Wood [58] also measured XPEO,water and obtained xpeo, water = 
1.35. Xu et al [55] measured the interaction between PS and PEO and reported Xps.peo = 0.02 ~ 
0.03. Zhu et al [5U] also measured Xps.peo an d and obtained Xps,peo = —0.00705 + 21.3/T. 
The polymerization index of PS-PEO diblock copolymer used in the experiments [5T] is typically 
~ 1000, thus we expect that xps,peo^ is sufficiently larger than the critical value XcN = 10.495 
[14] . Note that the polyelectrolytes such as poly(acrylic acid) (PAA) is often used as the hydrophilic 
subchian [49] [50] . Such polyelectrolytes can be dissolved into water easily, and thus we expect the 
X parameter between the polyelectrolytes and the water is negative. 

Thus we can say the diblock copolymers are in strong segregation region, and the effect of 
water addition is expected to be large for hydrophobic subchains, but not so large for hydrophilic 
subchains. Thus we consider that our simulation is qualitatively consistent with these experiments. 
Note that, however, we cannot compare the x parameter used in the simulations and measured by 
experiments directly because the value of the x parameter depends on the definition of the segment, 
and our free energy functional ([6]) is not quantitatively accurate as the SCF theory [22]. We also 
note that most experiments are carried out in the very strong segregation region in which it is 
quite hard to perform continuum field simulations. While it is difficult to compare our simulations 
with real experiments, our simulations are expected to give qualitatively correct physical process 
since our simulations take the crucial physical properties correctly; the hydrophobic interaction 
Xbs is sufficiently large and the hydrophilic interaction xas is small. These properties are not 
well considered in most of previous dynamics simulations. 

Other interesting properties are the size of the vesicles, the content of solvents inside vesicles, 
or the higher order structures. These properties cannot be discussed from our current results, since 
the system size is not large. We have only one vesicle in the simulation box (Figures 12131 and |4|) and 
there may be the finite size effect (see also Appendix [Bj . Nevertheless, we believe our model and 
simulation results are still valuable. These properties of vesicles should be observed if we perform 
the simulations for large systems while it is too difficult to perform large scale simulations by the 
current models, algorithms and CPU power. 

The characteristic time scales are also interesting and important property. As mentioned above, 
the real time scale can be estimated from the diffusion data. Here we estimate the time scale by 
using the empirical equation for diffusion coefficient of short polystyrene melts (at T g + 125°C) by 
Watanabe and Kotaka [BT] , 

D ~ 6.3 x 10- 5 M _1 cm 2 /s (27) 

where M is the molecular weight. 

We have to determine the number of monomers in one segment used in our simulations. Since 
the polymerization index of typical amphiphilic block copolymers used in experiments [50] [51] 
are about M ~ 50 ~ 1000, we consider that each segment used in simulations contains about 
10 monomer units. In our simulations, we used dimensionless segment size (b = 1). Thus we 
also need the size of the segment. The segment size of the styrene monomer is about 7A (it 
can be calculated from the experimental data of the radius of gyration [62] 63J). So if we assume 
the Gaussian statistics for monomers in the segment, size of the segment can be expressed as 
b ^ VlO x 7A. 
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Now we can estimate the characteristic time scale for the system. If we assume that the 
polymer chains are not entangled and obey the Rouse dynamics, the characteristic time scale r 
can be calculated as follows. 

_ 2C6g _ 2b 2 



k R T ND 



(28) 



60 is the size of the segment and N is the number of segments in one polymer chain (polymerization 
index), and D is the diffusion coefficient. N and M can be related as M = 104 x 10 x N and thus 
we have 

2 x 1040 x (VlO x 7 x 10- 10 m) 2 r , 

t~ ^ 5—5-; i-~1.6xl0 _6 s (29 

6.5xl0- 9 m 2 /s v ' 

Thus we know that the characteristic time scale in our simulations are about 1 microsecond, and 
the time scale of the vesicle formation is about 10 milliseconds. Of course the time scale estimated 
here is quite rough and cannot be compared with real experiments directly, but the time scale of 
our simulation is considered to be much smaller than one of experiments. The possible reasons are 
that the accuracy of the dynamic equation used for the simulation is not so good, and that the 
sizes of the formed vesicles are small since the system size is small. 



4.3 Comparison with the EPD simulations 

We compare our simulation results with EPD simulation results and discuss the reason why the 
EPD simulations [24j and our DF simulations give qualitatively different results (the mechanism 
II and the mechanism I). 

The main difference between the EPD simulations and our DF simulations are the free energy 
functional model and the dynamic equation. The EPD simulations use the free energy calculated 
from the SCF theory while the DF simulations use the free energy functional ©■ It is known 
that the DF theory is less accurate than the SCF theory, but gives qualitatively acceptable results 
[22l [18] . Thus it is difficult to consider that the accuracy of the free energy model affect the 
observed mechanisms. We expect that the mechanism I is reproduced if we change the free energy 
functional model from eq © to more accurate SCF model. 

Other difference is the dynamic equation. The EPD model employs nonlocal mobility model, 
but in the derivation of the dynamic equation, rather rough approximations are involved 64] . 
Because of these approximations, the local mass conservation is not satisfied in the EPD simulations 
(unless the system is homogeneous). This may affect the dynamic behaviour seriously. The mobility 
in our dynamic equation ([8} is the local type, and thus it is considered to be less accurate than the 
nonlocal model. But in our dynamic equation, local mass conservation is satisfied exactly. Thus 
we consider that the EPD simulations reproduce the mechanism II because it does not satisfy the 
local mass conservation. We expect that if the dynamic equation which satisfy the local mass 
conservation is used with nonlocal mobility and the SCF free energy, the mechanism I should be 
reproduced. 



5 Conclusion 

We have shown that we can reproduce the vesicle formation mechanism I by the simulation using 
the dynamic density functional equation (|8|) and the free energy functional |6]). The simulation 
results are qualitatively in agreement with other simulations [3l [4] . This is the hrst realistic vesicle 
formation dynamics simulation based on the field theoretical model. To reproduce the mechanism 
I, we need the free energy functional which can be applied to strong segregation regions, the large 
thermal noise, and the proper interaction parameters for hydrophilic and hydrophobic interactions. 
These conditions are not satisfied previous field theoretical dynamic simulations. We have also 
shown that the formation dynamics of various micellar structures (spherical micelles, cylindrical 
micelles, vesicles, bilayers) and morphological transitions can be simulated by our model. The 
morphological transitions can be reproduced only by changing the x parameter which corresponds 
to the change of the solvent condition. 
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The merit of using the field theoretical approach is that we can use the Flory-Huggins interaction 
parameters, instead of potentials between segments which is required for particle simulations. 
We have shown that we can mimic the solvent condition control by changing the \ parameter 
between the hydrophobic segment and the solvent and have shown that the solvent condition control 
actually causes morphological transitions. For coarse grained particle simulations, we cannot use 
the microscopic potential directly, and in most cases the phenomenological potential models are 
used. Of course for some particle simulations, such as the DPD, we can use the \ parameters by 
mapping them onto the DPD potential parameters. But such method is based on the parameter 
fitting [65] and not always justified. Thus we consider that especially for many component systems, 
such as solutions of amphiphiles, the field theoretical approaches have advantages. 

Although the x parameters in the current model cannot be compared with the experiments 
directly, our data may help understanding experimental data or the physical mechanism. It will 
be possible to use the experimentally determined x parameters by employing more accurate free 
energy model such as the SCF model [151 EE HZ] and efficient numerical algorithms. (We note 
that recently Honda and Kawakatsu [66] proposed the hybrid theory of the DF and the SCF 
and reported that accurate and numerically efficient simulations can be performed by the hybrid 
method. By employing their method, accurate and fast simulations for the dynamics of micelles 
and vesicles may be possible.) The quantitative simulations will be the future work. 

While the efficiency of our simulation method based on the density functional theory (eqs ([6]) 
and (fl4|) ) is larger than the SCF simulations, currently it is still comparable to particle simulations. 
Besides, because the late stage of the mechanism I is a very slow process, the current model still 
need considerable computational costs. To study the late stage dynamics or large scale systems (for 
example, many vesicle systems), we will need more coarse grained and numerically more efficient 
computational methods. Our current work can be the basis of the further coarse grained field 
theoretical model. 
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Appendix 

A Noise Generation Scheme 

In this section we derive the numerical scheme for generating random noise field which satisfies the 
fluctuation dissipation relation (eqs (jTSj) and (p~6|) V Such a noise field can be expressed as [42] 




(30) 



where Ui(r,t) is the Gaussian white noise (vector) field which satisfies 



(«i(r,t)> = 
( Wi (r, QwjCr'.F)) = 5 l3 5{r - r')S(t-t')l 



(31) 
(32) 



where 1 is the unit tensor. 
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To generate the noise numerically, we need the descretized version of eq ([50)1 . The descretized 
noise field £i(r,i) is defined only on the lattice points. Here we express the position of the lattice 
point as r — n x h x + n y h y + n z h z where h a (a = x, y, z) is the lattice vector and n a is the integer. 
The rescaled time t is also descretized as t = n t At where n t is the integer. The Laplacian operator 
is replaced by the standard centre finite difference operator as follows 



V 2 /(r)^ T^lf(r + h a )~f(r) + f(r-h a )] 



(33) 



where f(r) is a descretized field defined on the lattice points. The delta functions are replaced as 
follows 



S(r-r') 



\h x \\h y \\h z 



At 



(34) 
(35) 



where 5 rr > = ^.ti^n^^n^; and %, = 5 n «nj- 

The fluctuation dissipation relation for the second order moment (eq (11611 s ) can be rewritten as 



Zi(rM(r',i')) = • [tf(r)V6(r - r')] S(t-t') 



, 2 , ^ipi{r)S(r - r') 



#(r)V- 



s(t-t') 



-4/3-%- [^(r)V 2 [^(r)5(r - r')] - bPi(r)S(r - r')] V 2 ^(r)] S(t - f ) 



(36) 



Using eqs (l33l) - (|35|) . we get the descretized version of the fluctuation dissipation relation ([36 



^i(r)^j(r - h Q ) [5 Pr / - <J( r -h )r'] 



(37) 



In this work, we employ the following descretized equation to generate the noise field. 

1 



a—x,y,z 



(38) 



y/Mr)Mr - h a ) 4 a \r - h a /2,i) 



where uj^ a \r,t) is the a element of U3i{r,t). w,-(r,t) is the descretized Gaussian white noise field 
which satisfies 



(wi(r,t)) =0 



$i'i$rv f &+¥> 



-a 



/iJ|/iJ|hJAt' 



(39) 
(40) 



Notice that ui,-(r , ,t) is defined on the staggered lattice. The position of the staggered lattice point 
is r = (n x + l/2)h K + n y ^+n z ^ z ,n x h x + (n a + l/2)h I/ + n z ft. z ,n x /i x +n a hj / + (n z + l/2)h z where 
n a is the integer. It is easy to show that the noise field generated by eq ([38| satisfies eq ([37[) . 
The noise generation scheme is finally written as follows. 

1. Generate the Gaussian white normal distribution vector noise field u>i(r,t) on the staggered 
lattice points. The noise is generated by using the Mersenne twister pseudo-random number 
generator [57] and the standard Box-Muller method [4"rj] . 



2. Calculate the noise field £i(r,t) from ipi(r) and u>i(r, t) by using eq ([38 
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B Finite Size Effect 



It is well known that the size of the simulation box often affects the simulation results (the finite 
size effect). In this appendix, we study the finite size effect for our simulations by performing 
simulations with different box size. 

Figure |TD] shows the simulation results for different (small) box sizes. All the parameters except 
for the box size are the same for the previous vesicle formation simulation (Figure [2]) . The box 
sizes are set to 66 x 66 x 66 (12 x 12 x 12 lattice points), 86 x 86 x 86 (16 x 16 x 16 lattice points), 
126 x 126 x 126 (24 x 24 x 24 lattice points), and 166 x 166 x 166 (32 x 32 x 32 lattice points). 
We can observe the spherical micelles (Figure [TOT a)) and the cylindrical micelles (Figure [TDT b) .(c)) 
for small boxes. These structures are much different from the vesicle structures in Figure [TOTd) or 
Figure [5] This means that if the simulation box is too small, we cannot obtain vesicles. 

It also implies that there may be the finite size effect in our vesicle formation simulations (box 
size 48 x 48 x 48), since the size of obtained vesicles are in most cases comparable to the box sizes. 
But at least we can claim that our vesicle formation process itself is not affected qualitatively 
since the system size is not so small for the formation of disk like structures (unlike the case of 
Figure [Tola) . (b) . (c) ) and we can observe that the small vesicle formation process (Figure E^d)) is 
qualitatively same as one for larger vesicles. 

C Cutoff Length Dependency 

The cutoff length A for the long range interaction term in eq is introduced rather intuitively, 
and therefore the value of A for simulations should be considered carefully. In this section, we 
perform the simulations with different values of A. Simulations are performed in two dimensions. 
Parameters are as follows: N = 10, f A = 1/3, /b = 2/3 A = 2,4,5,6,8. 4> s = 1 - 4> p = 0.8. 
Xab = 2.5,% J 45 — — 0.5,xbs = 5, b = 1. The size of the 2 dimensional simulation box is 
646 x 646 (128 x 128 lattice points, periodic boundary condition). The time step is At = 0.0025 
(for A = 2, 4, 5), 0.000625 (for A = 6), At = 0.0003125 (for A = 8). The magnitude of the noise is 
/3 _1 = 0.078125. All the simulations are started from the homogeneous state. 

Figure [Til shows the results for various value of the cutoff length ((a),(b) A = 1, (c),(d) A = 2, 
(e),(f) A = 4, (g),(h) A = 5, (i),(j) A = 6, and (k),(l) A = 8) at t = 312.5, 1250. We can observe the 
droplet like structures are observed for A = 1. and micellar structures for 2 < A < 8. The resulting 
morphologies for 2 < A < 8 are qualitatively same. (Strictly speaking, A = 2 looks slightly different 
from others and somehow resembles to A = 1.) As expected, the cutoff length dependence of the 
morphology is small. Thus even if the value of A is not accurate, the results are expected not to be 
affected qualitatively. From this result, we can justify to use the value A = 5, which is estimated 
by the rough argument. 
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Figures 



Figure Captions 

Figure [JJ Schematic representation of two vesicle formation mechanisms from the initial ho- 
mogeneous state. Black and grey color corresponds to hydrophobic and hydrophilic subchains, 
respectively, (a) mechanism I: First, small micellar structures are formed. The micellar structures 
grow by collision and become large cylindrical or open disk-like micelles. Finally the open disk-like 
micelle close up to form closed vesicles, (b) mechanism II: The initial stage is similar to the mech- 
anism I. However, the small micellar structures formed in the initial stage grows up to be large 
spherical micelles. The large spherical micelles then take the solvent inside and thus vesicles are 
formed. 

Figured) Snapshots of the dynamics simulation for the amphiphilic diblock copolymer solution 
(<f>p = 0.2). Grey surface is the isodensity surface for the hydrophobic subchain, 0s(r) = 0.5. (a) 
t = 6.25, (b) t = 312.5, (c) t = 1562.5, (d) t = 3125, (e) t = 4687.5, (f) i = 6250. 

Figure [3] Snapshots of the dynamics simulation for the amphiphilic diblock copolymer solution 
((j) p = 0.2). Grey surface is the isodensity surface for the hydrophobic subchain, 0b (r) = 0.5. (a) 
t = 6.25, (b) t = 312.5, (c) t = 1562.5, (d) t = 3125. 

Figure |U Snapshots of the dynamics simulation for the amphiphilic diblock copolymer solution 
((j)p = 0.2). Grey surface is the isodensity surface for the hydrophobic subchain, 0s(r) = 0.5. (a) 
t = 6.25, (b) t = 3125, (c) t = 4687.5, (d) t = 6250. 

Figure [5j Two dimensional cross sections of density field for (a) 4>A{f) and (b) ^s(r). (c) The 
density correlation function Sab(q)- t = 6250. 

Figure [6) Snapshots of the dynamics simulations for amphiphilic diblock copolymer solutions 
with various volume fractions. Grey surface is the isodensity surface for the hydrophobic subchain, 
4> B {r) = 0.5 at t = 3125. (a) 4> p = 0.1, (b) 4> p = 0.15, (c) 4> p = 0.25, (d) 4> p = 0.3. (</> p = 0.2 
corresponds to Figure [^d).) 

Figure [7j Snapshots of the dynamics simulations for amphiphilic diblock copolymer solutions 
with various xbs (the \ parameter between the hydrophobic chain and the solvent) . Grey surface 
is the isodensity surface for the hydrophobic subchain, 0b (r) = 0.5 at t — 3125. (a) xbs — 2, (b) 
Xbs = 2.5, (c) xbs = 3, (d) xbs = 3.5. (xbs = 5 corresponds to Figure H^d).) 

Figure [H] Snapshots of the structural transition dynamics simulations. The x parameter xbs 
is initially set to xbs — 5 and at t — 6250, xbs is set to lower value. Grey surface is the 
isodensity surface for the hydrophobic subchain, 0b (r) = 0.5 at t = 9375. (a) xbs = 5 — > 2.5, (b) 
Xbs = 5 — ► 3. 

Figure [9] Snapshots of the dynamics simulations for values of /3 _1 at t = 6.25,12.5. (a),(b) 
/H = 0.00390625, t = 6.25, 12.5. (c),(d) /3" 1 = 0.0390625, i = 6.25, 12.5. (e),(f) = 0.390625, 
t = 6.25, 12.5. 

Figure 1101 Snapshots of the dynamics simulations for amphiphilic diblock copolymer solutions 
with various small box sizes at t = 3125. (a) box size 66 x 66 x 66 (12 x 12 x 12 lattice points), 
(b) box size 86 x 86 x 86 (16 x 16 x 16 lattice points), (c) box size 126 x 126 x 126 (24 x 24 x 24 
lattice points), and (d) box size 166 x 166 x 166 (32 x 32 x 32 lattice points). Super cells (size 
246 x 246 x 246) is shown. Small gray boxes show the real simulation box. 

Figure [Til Snapshots of the dynamics simulations for amphiphilic diblock copolymer solutions 
with various cutoff length A. Black color represents 0b(»")- (a)j(b) A = 1, t = 312.5, 1250. (c),(d) 
A = 2, t = 312.5, 1250. (e),(f) A = 4, I = 312.5, 1250. (g),(h) A = 5, t = 312.5, 1250. (i),(j) A = 6, 
t = 312.5, 1250. (k),(l) A = 8, t = 312.5, 1250. 
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(a) Mechanism I 



(b) Mechanism II 

Figure 1: 
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Figure 6: 
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Figure 9: 
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Figure 11: 
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